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xl ! Abstract 

Space and time scales are not independent in diffusion. In fact, numerical simulations 
^ . show that different patterns are obtained when space and time steps ( Aa; and At) are varied 

independently. On the other hand, anisotropy effects due to the symmetries of the dis- 



O . cretization lattice prevent the quantitative calibration of models. We introduce a new class 

of explicit difference methods for numerical integration of diffusion and reaction-diffusion 



equations, where the dependence on space and time scales occurs naturally. Numerical 
solutions approach the exact solution of the continuous diffusion equation for finite Ax 
c/2 I and At, if the parameter 77V = -DAt/(Ax)^ assumes a fixed constant value, where A^ is an 

odd positive integer parametrizing the alghorithm. The error between the solutions of the 
discrete and the continuous equations goes to zero as (Ax)^*-'^"'"^-' and the values of 7_;v are 
dimension independent. With these new integration methods, anisotropy effects resulting 
^ ', from the finite differences are minimized, defining a standard for validation and calibration 

of numerical solutions of diffusion and reaction-diffusion equations. Comparison between 
numerical and analytical solutions of reaction-diffusion equations give global discretization 
errors of the order of 10~^ in the sup norm. Circular patterns of travelling waves have a 
maximum relative random deviation from the spherical symmetry of the order of 0.2%, 
and the standard deviation of the fiuctuations around the mean circular wave front is of 
the order of 10""^. 
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1 - Introduction 

After the important work of Turing [1952] on ttie ctiemical basis of morphogenesis, 
theoretical and experimental studies have shown that solutions of reaction-diffusion (RD) 
partial differential equations present self-organizing properties, common to biological and 
chemical extended systems. These properties appear as emerging coherent patterns or 
structures in spatially extended media. Due to the nonlinear nature of local kinetic me- 
chanisms in RD systems, local fluctuations can grow and propagate by diffusion to the 
surrounding media. For this reason, extended media supporting local nonlinear kinetic 
mechanisms are called excitable media [Winfree, 1990; Keener & Tyson, 1992]. 

The prototype experimental system for the study of pattern formation in extended 
media are the Belousov-Zhabotinsky reaction [Zaikin & Zhabotinsky, 1970; Zhabotinsky 
& Zaikin, 1973] and, in biological systems, the aggregation patterns of colonies of Dic- 
tyostelium discoideum, [Tyson et al., 1989; Steinbock et al., 1993]. Due to the similarity 
of dynamic behaviors, these experimental systems provide insights into the mechanisms 
of pattern formation (morphogenesis). A large collection of experimental evidence of spa- 
tiotemporal patterns in nature can be found in Field, Koros & Noyes [1972], Meinhardt 
[1982], Harrison [1993], Cross & Hohenberg [1993], Kock & Meinhardt [1994] and Epstein 
& Showalter [1996]. 

Reaction-diffusion (partial differential) equations are the simplest mathematical mod- 
els for the study of pattern formation in excitable media. The coherent patterns appearing 
in numerical simulations of RD equations are undamped travelling waves, spiral travelling 
waves, stable strips and spotty structures (Turing patterns), [Castets et a/., 1990]. From 
the qualitative point of view, the same structures can be obtained with cellular automata 
models, [Markus & Hess, 1990; Durrett & Griffeath, 1993]. However, cellular automata 
models cannot be calibrated and validated from real parameters (kinetic reaction rates 
and diffusion coefficients), as the physical and chemical mechanisms that characterize the 
observed phenomena are generally continuous models. These parameters are specially 
important as different pattern topologies are obtained when they are varied (bifurcation 
phenomena) . 

The possibility to validate a physical or chemical mechanism for the emergence of a 
given coherent pattern lies on the existence of reliable numerical alghorithms. If simula- 
tion models show the same qualitative behavior as real systems, model parameters can be 
calibrated with experimental data. If the calibrated model has the same behavior as the 
real system, the model is validated, and the association of specific parameter values to the 
real system characterizes its properties, [Oreskes et al., 1994]. A classical example of this 
methodology is obtained in the study of diffusion. In this case, the diffusive properties of 
a medium are quantified by the diffusion coefficient, measured by comparing experimental 
concentration profiles with the analytical solutions of the diffusion equation. Therefore, 



diffusion, originated by the random impacts of the solvent molecules with suspended par- 
ticles, is quantifiable by the diffusion coefficient, measuring the mean square displacement 
of the suspended particles [Chandrasekhar, 1943]. 

As experimental RD systems are in general nonlinear, the calibration and validation 
of models should be based on the comparison between the data obtained with numerical 
solutions of RD equations and data acquired in experiments. This rises the problem of 
comparison between numerical and analytical solutions of nonlinear equations — bench- 
marking. 

1.1 - Problems with numerically obtained patterns with RD equations 

The prototype mathematical model for the study of pattern formation in excitable 
media is the two dimensional system of RD equations 

where (p = ((/?i,(/?2) are the concentrations of some chemical species or morphogens [Tur- 
ing, 1952], D^ = (1)1,1)2) are diffusion coefficients and X'^ = (Xi,X2) is a vector field 
describing the local kinetics of the system. Systems of partial differential equations of this 
type describe, for example, the space-time pattern formation in the Belousov-Zhabotinsky 
reaction and the interaction between cells separated by permeable membranes, [Turing, 
1952; Zhabotinsky & Zaikin, 1973; Murray, 1993]. 

In general, for non-linear vector fields, there are no general analytical techniques to 
obtain solutions of (1.1) and we must rely on numerical methods. Stability and consistency 
conditions for finite difference integration methods insure that the truncation errors of 
numerical solutions vanish when the discrete space and time steps go simultaneously to 
zero [Smith, 1985; Sewell, 1988]. However, for large integration times the accumulated 
error can grow indefinitely [John, 1971], distorting concentration patterns. 

The simplest finite difference explicit method used in numerical integration of eq. (1.1) 
is the system of difference equations 

if\+f = AtXi(vpl,,,,, ifil,,^) + -^^ (y'U-i.j + y'U+i,, + y'U,,-! + v^U.^+i - 4 ^Uj) 

(1.2) 
where the subscripts (z,j) refer to the coordinates (iAx^jAx) of the vertices of a squared 
symmetric lattice in the plane. The solutions at time t, y^i ^^j and ^2,i,j^ ^^^ obtained 
iteratively through given initial functions (/7i(x, y, 0) = /i(x, y) and ip2{.x^y,Q) = f2{x,y), 
Vi,i,j = /i(^^3:,jAx) and </'2,i,j — f^iiAx^jAx). The centered difference scheme (1.2) is 



stable if 7 = maxjDi, D2}At/(Ax)^ < 1/4, and the spatial truncation error is of the order 
of (Ax)2, [Smith, 1985; SeweU, 1988]. 

In the literature on numerical methods, it is generally accepted that for the integration 
of two-dimensional diffusion and RD equations, the parameter 7 = D /S.t / {/S.x)'^ must be as 
close as possible to the value 7 = 1/4, [Carslaw & Jaeger, 1959; Crank, 1975; Smith, 1985; 
SeweU, 1988]. This is justified by the simple fact that computation time is faster and the 
choice of 7 is not very important if the purpose is to determine equilibrium solutions [Press 
et al, 1989]. In fact, if X = 0, equilibrium solutions of the difference scheme (1.2) coincide 
with the equilibrium solution of (1.1). However, if the goal is to obtain information about 
propagation velocities of diffusion fronts, or to study the pattern formation in reaction- 
diffusion processes, this approach presents serious problems, as we show now. 

The first observation is that numerically computed concentration profiles obtained 
with RD equations depend strongly on the magnitude of 7 = DAt/{Ax)'^, within the 
stability condition 7 < 1/4. To show this, we take the Brusselator model, [Prigogine & 
Lefevre, 1968], in an extended two dimensional media, defined by the equations 

-df ^ ^^''^^ ' ^^^^^^ + '"'[-8^ + -W) 

where the kinetic terms correspond to the system of autocatalitic reactions A -^^^ t/Ji, 
B + ifi -^^'^ {p2 + -D, 2(/7i + {p2 -^^^ 3(/?i and ipi -^^'^ E, the kiS are the reaction rates, and 
the concentrations A and B are kept constant. System (1.3) has an equilibrium unstable 
solution for ipi = Aki/k^ and (/?2 = Bk2k4/{Akiks). 

In Fig. 1, different evolved patterns are shown, obtained with numerical integration 
method (1.2), for the same initial condition on a lattice of 200 x 200 cells and several values 
of 7 = max{Di, D2}At/(Aa;)^. The initial condition has been chosen to produce spiral 
patterns, [Klevecz et a/., 1992]. As 7 increases within the stability condition 7 < 1/4, pat- 
terns change shape and, after some threshold, the spiral disappears leading to concentric 
waves. For small values of 7, the actualization process induces a squared spurious sym- 
metry in the concentration profiles, showing spatial anisotropy. These simulations show 
that different nonequilibrium patterns are obtained without changing either the physical 
parameters (diffusion and kinetic coefficients), or the initial conditions. On the other hand, 
as Ax = a/ -D At/7, different lattice lengths are obtained as we vary 7, for the same number 
of lattice sites. 

Therefore, the propagation velocity and physical dimensions of patterns obtained nu- 
merically, as well as its symmetry properties, are strongly dependent on the magnitude 
of 7, leading to deformations in the patterns, together with quantitative changes in the 
propagation velocities of wave fronts. 



On the other hand, the ad hoc scahng introduced in Fig. 1, show that space and 
time variables are not independent. This is a well known symmetry associated with the 
diffusion equation, which is invariant under the linear substitutions x' = ax and t' = a^t 
[John, 1971]. 

So, we are faced with the problem of knowing which of the simulations of Fig. 1, 
better approximates the solution of the continuous RD equation, for the prescribed initial 
data. 

To avoid the problem of numerical diffusion anisotropy in patterns of RD systems, 
some authors use implicit integration schemes as the Crank-Nicholson or ADI methods, 
[Press et al., 1989], which are computer time consuming. With explicit integration discrete 
methods, Barkley [1990] proposed an explicit nine point formula difference approximation 
for the Laplace operator, Pearson [1995] adds random noise, and Weimar et al. [1992] 
analyse several "masks" for the discrete approximation to the diffusion equation. In the 
context of ecological modelling, Dejak et al. [1987] makes an error analysis as a function 
of the scaling parameter 7. For cellular automata models, Markus & Hesse [1990] use a 
random algorithm. All these strategies to overcome the problem of diffusion anisotropy are 
heuristic and estimates of global numerical errors are lacking. On the other hand, these 
methods do not show the role of the heuristic scaling relation introduced in depicting the 
simulations of Fig. 1, which otherwise change the physical dimensions of the system. For 
a general discussion on "plausible-looking results which are qualitatively incorrect" see 
Ruuth [1995]. 

Another problem when simulating diffusion processes with explicit methods is related 
to the maximum propagation velocity of "suspended particles" . For the discrete method 
(1.2), in one time step, local concentration changes are due to next neighbors contributions 
and the maximum propagation velocity is Vmax = 2^. But, for the diffusion equation, 
Vmax = +00, as is Well know. The solution of the discrete explicit method (1.2) converges 
to the solution of Eq. (1.1), if limAt^o hiHAx^o 7 =constant [John, 1971]. But as 7 = 
-DAt/(Ax)^, this last condition implies that limAt^o limAx^o ^max = +00, but for finite 
Ax and At, Vmax remains finite. Therefore, the best choice of 7 remained an open question. 

According to the microscopic interpretation of diffusion by the theory of Brownian 
motion, the value of 7 should be, 7 = 1/2, 7 = 1/4 and 7 = 1/6 in one, two and three 
dimensions, respectively [Murray, 1993]. So, it seems natural to choose one of the above 
values of 7 in simulations. However, these values correspond to the upper limit of stability 
of explicit methods in one, two and three dimensions, respectively, and introduce numerical 
instabilities due to roundoff effects. 



1.2 - An overview of results and organization of the paper 

In this paper, we derive a new class of explicit difference methods for the numerical 
integration of diffusion and reaction-diffusion equations and a rigorous criterion for the 
choice of the space-time scaling parameter 7 = D /S.t / {/S.x)'^ in order that: 

a) Explicit difference methods for the numerical integration of RD equations give accu- 
rate results for the transient solutions of diffusion and reaction-diffusion processes, 
enabling to calibrate model parameters with experimental data. 

b) The relation between space and time scales, limAt^o liniAx^o 7 =constant, a property 
of diffusion equation, is correctly interpreted. 

c) The spatial distribution of concentrations do not depend on the symmetries of the 
discretization lattice. In particular, diffusion anisotropy in dimensions two and three 
are minimized, without additional hypotheses or data modification. 

The calibration of the space and time scales results from a choice of a constant value 

for 7, dependent on the order of the global discretization error: 7 = 7Ar, for an error of 

the order of {/S.x)'^^^^'^\ where A^ > 1 is an odd integer. 

The main consequences of the approach developed here are: 
i) There exists an optimum value of 7, 7 = 77V, for which the new class of explicit 
difference methods shows the smallest global error, when its solutions are compared 
with the solutions of the continuous diffusion and RD equation — Fig. 2. The value 
of 7Ar is dimension independent, defining a space-time scaling relation, 
ii) The algorithm is explicit and the accuracy of the simulations can be arbitrarily in- 
creased, for finite Ax and At — Eqs. (2.6), (2.34) and (2.47). 

iii) Numerical lattice anisotropy in diffusion and reaction-diffusion equations are mini- 
mized — Fig. 5. 

iv) In the context of the probabilistic interpretation of diffusion through Brownian mo- 
tion, the transition probabilities from one lattice cell to the adjacent cells is exactly 
calculable, being dependent on the number of neighborhood cells taken into account in 
the discretization. The values of the transition probabilities determine the optimum 
Laplacian mask operator [Weimar et al., 1992], as a function of neighborhood order. 
v) Solutions of RD equations with cylindrical or spherical symmetries numerically coin- 
cide, independently of the space dimension of simulations — Fig. 8. 

In the next section we introduce the main technique for the discretization of the 
diffusion equation in dimensions 1, 2 and 3, in order to account for lattice anisotropies. 
The main results of the paper are stated in Theorems A, B and C, and the exact solutions 
of a linear RD equation are compared with the numerically calculated solution — bench- 
marking. In Sec. 3 we summarize the main results of the paper from the applied point of 
view. 



2 - Main Results 

2.1 - One space dimension 

We begin with the diffusion equation in one space dimension 

^-D^ (2 1) 

dt~^dx^ ^^-^^ 

and we take a lattice E in the {x,t) slab with x G R and < t < T. With Ax and 
At as the increments in the variables x and t, the lattice E has coordinates {nAx,kAt), 
with n G Z and /c G N. If (j){x,t) is a solution of the parabolic equation (2.1), satisfying 
the initial condition (j){x,0) = /(x), this solution, restricted to E, will be denoted by 
(f)^ = (t){nAx,kAt). 

As it is well known, for the initial condition /(x), the solution of the diffusion equation 
(2.1) is 

Hx,t)= f(y)^=^e-^'ii^dy (2.2) 



If f{x) is a discontinuous and integrable, the function (j){x,t) is of class C°°(R), for every 
t > 0. So, without loss of generality, we consider that (/){x,t) is of class C°°(R). On the 
other hand, if f{x) has compact support, 4>{x^t) ^ for every t > and finite x. This 
implies that diffusive propagation has infinite velocity. 

In order to discretize the diffusion equation (2.1) on the lattice E, we introduce the 
usual finite difference approximation to the space and time derivatives, 

d(j) (j){x,t + At)-(/){x,t) 



dt At 

d'^(j) (p{x — iAx, t) + (j){x + zAx, t) — 2(p{x, t) 



2. ... .A..^ , ... , .A..^ o... .^ (2-3) 



dx'^ i'^{Ax)'^ 

where i is some positive integer. To account for large propagation velocities in the dis- 
cretization of the diffusion equation, we rewrite (2.1) in the form 



^ Q2, 



i-^E«.§J (2.4) 



dt ^^-^ dx 

1=1 



where the a^ are non negative constants such that 



TV 

5^a^ = l (2.5) 

1=1 



Introducing the finite difference approximation (2.3) into (2.4), we obtain the discrete 
equation 



N 
V. 

i=l 



r^ - ^n + 7iV E S i<-^ + <+r - 2^') (2-6) 



where 

'" = °(A^ ("> 

The integer N is the space width of the finite difference equation (2.6) and the constants 
ai measure the connectivity strength between each ceU and its neighbors, up to width A^. 
With v^ = (j)'^ = f{nAx)j aU the v'^, for /c > 0, are determined recursively. 

For A^ = 1, the relation between the recursive solutions of (2.6) and the solutions 
of the diffusion equation for the same initial function f{x) = (p{x,0) is easily obtained 
through the maximum principle. As a matter of fact, if 7Ar is a sufficiently small constant, 
then 

llk(^,i)lll< 11/(^)11 (2.8) 

where 

\\f{x)\\ = supx\f{x)\ and \\\v{x,t)\\\ = supx^t\v{x,t)\ 

is the sup norm taken for x G R and < t < T. Under these conditions, if Ax, At -^ 0, 
v'^ converges uniformly to the solution of the diffusion equation. In the case A^ = 1, the 
uniform convergence is achieved under the stability condition 71 = DAt/{Ax)'^ < 1/2, 
(Lemma II, §7.2 of [John, 1971]). 

For the general case of A^ > 1, we apply Schwartz inequality to (2.6): 

N N 

i=l i=l 

So, if < 7Ar X]n=i ^ — ^ ^^^ ^^^ constants a^ are non negative, we have 

\\Un II ^ llt^nll l^-^J 

With \\v':^\\ < ||/(a^)|| and iterating (2.9), we obtain inequality (2.8). Therefore, the finite 
difference method (2.6) is stable and obeys the maximum principle if 

and ai > 0, for i = 1,...,N. By the Lax Equivalence Theorem, stability of (2.6) is a 
necessary and sufficient condition for the initial-value problem of the diffusion equation to 
be properly posed [Richtmeyer & Morton, 1967]. 

However, from the numerical point of view there is no control on the discrete incre- 
ments Ax and At, or on the values of the parameter 77V- The next theorems solve these 
problems showing that convergence of the solutions of finite difference diffusion equation 
to the solutions of the diffusion equation is obtained in the limit N ^ 00, for finite Ax 
and At. To be more precise, we denote by v^ [x^t] the solution of the difference equation 
(2.6) on the lattice E, for some choice of the width A^. Our main objective is to establish 
the conditions under which limAr^oo v^ {x, t) = ^(x, t), for finite Ax and At. 



Theorem A. If N is an odd positive integer, then there exists a positive constant 'Jn, 
and real constants cti, . . . , a^, solutions of the system of equations 



Z! 



(2i)\ ^ 



I'i 2 



AT 



— yj]«.n2^ = (2.116) 



(iV + 1)! (2(iV + l))!^^ 

Moreover, if ai, . . . ,aN are non negative and '-jn < 7jV) with •^'^ given by (2.10), the 
solution of the difference equation (2.6) approaches the solution of the continuous equation 
(2.1) with an error of the order of (Ax)'^^^'^'^' . The values of ai as a function of "fN are 
obtained solving the system of linear inhomogeneous equations (2.11a) and 'Jn is largest 
real root of the polynomial (2.11b). 

Remarks on Theorem A: As a matter of fact, we have tested the positivity of the 
constants ai, as weU as the stability condition 77V < 7jv5 up to A^ = 23, for odd A^. 
Numerical tests indicate that the constants ai are positive for A^ odd, only if •^n is the 
largest real root of (2.11b). This suggests the stronger assertion that, lim^^oo v'^^~^^{x, t) = 
4>{x, t). For A^ even the theorem is false. 

Proof: Let (j){x, t) be a solution of the diffusion equation (2.1). We suppose that (/){x, t) 
is of class C°°(R) in both variables x and t. We define the difference operator 



N 

a 



A(j) = (j){x, t + At)- (j){x, ^) - 77V X] ^ {(f){x + nAx, t) + (^(x - nAx, t) - 2(j){x, t)) (2.12) 



n 

where X]n=i ^n = ^ and 7Ar = DAt/{AxY. Applying the Taylor theorem to the operator 
A(/), we obtain, 

i=l \ ^ ' n=l 

N 



^ dx^i^+^)^ ' ^""{{N + Dl {2{N + l))\^^^'' ^+uy^l^x) 

(2.13) 



n=l 



where we have used the relation, 



d'-cj) j d'^\ 



Choosing 77V and «!,...,«„ such that 



^-7^J2^r.n'(^~'^=0, t = l,...,N (2.14a) 

il (2zj! ^ — ^ 



n=l 



(^-(2(^E-""-0 (2.1«) 



n=l 



we have, for large A^ and sufficient small Ax, |A(^| = O ((Ax)^*-'^"'"^-'). 

Let v^ be a solution of the finite difference equation (2.6), and let A^v be the operator 
defined by 



N 

ar 



A^t^ = vt^' -v^-1nJ2^ (^'+^ + 'In - ^vf) (2.15) 

As, A^v = 0, we have 






\\AcP - A^t;|| = O ((Aa:)2(^+2)) (2.16) 

With, ef = (p{iAx, kAt) — v^ and subtracting (2.15) from (2.12), 



N 



A<t> - A^v = er ^ - ef - 7N E S (^'+n + el^ - 2ef ) 
and 






N 



e1+' = A<p- A^v + e,^ + 7^ E S (e'+n + et. - 2e,^) 
Therefore, by the Schwartz inequality, 

\\e'^+^\\ = \\(j){x,t + At)~v^ix,t + At)\\< \\A(l)-A^v\\ + mx,t)-v^{x,t)\\ (2.17) 

on the lattice E, provided a^ > 0, for z = 1, . . . , A^, and < 7jv < 7jv- As, ||(/)(x, 0) — 
t;^(x, 0)11 = in E, and iterating (2.17), 

\Mx,t) - z;^(x,t)||| < ^^O ((Aa:)2(^+2)) = O ((Aa:)2(^+2)) (2.18) 

for all (x, t) G E. 

To finish the proof, we investigate the conditions under which equations (2.14) have 
real and positive solutions. The system of A^ equations (2.14a) can be written in the form 

Aa = L (2.19) 

where A is a Vandermond matrix with a^j = X^~ = (j^)*~^, and the vectors a and L have 
components ct^ and (2i)!7]y"^/2(i!), respectively. So, as det A = ni<7(-^j ~ -^0 ¥" 0? the 
equation (2.19) has solutions a^ = ai(7Ar), for every A^ > 0, where ai(7Ar) are polynomials 
in 7jv. Introducing the polynomials ai(7Ar) into (2.14b), and as the ai(7Ar) have degree 
A^ — 1 in 7jv, (2.14b) has a real root if A^ is odd. 

We now prove that (2.14b) has a positive real root. With (2.14a) for i = 2, . . . , N and 
(2.14b), we construct the linear equation 

Bd = L' (2.20) 

10 



where S is a Vandermond matrix with elements bij = A* = (j^)\ and L' has components 
(2z + 2)!7]y/2((z + 1)!). As detS 7^ 0, the solutions of equation (2.20) are polynomials in 
■jn of degree A^. We denote these polynomial by a^(7Ar). Introducing these solutions into 
(2.14a), for z = 1, we have 

N 

p{^^):=J2<hn)-l = (2.21) 

i=l 

As detS 7^ 0, solving (2.20) for 7^ = 0, we have a[{0) = 0, and p{0) = —1. Therefore, if 
A^ is odd, p{'^n) has a positive real root. Thus, Theorem A is proved. 

In Tab. 1 we present the values of 7Ar and a^, for A^ = 1,3 and 5, determined by 
Theorem A, as well as the value of the parameter 7^, defining the stability limit of the 
difference equation (2.6). 

Theorem A has several important consequences for the numerical simulation of diffu- 
sion equations, namely: 

1) The constant 7Ar = DAt/{Ax)'^ introduces a coupling between the space and time 
scales, showing that in the finite differences approximation to the diffusion equation the 
space and time scales are not independent, for A^ odd and finite Ax and At. As 7Ar is 
a constant, it defines a scaling relation for the discrete diffusion processes. This scaling 
relation is the best one in the sense of error minimization in the sup norm. 

2) The global error between the numerical and exact solutions of the diffusion equa- 
tions obtained by the finite difference method (2.6) depend only on the space increment Ax 
and is time independent. Therefore, choosing A^, the accuracy of the numerical solutions 
can be increased by decreasing Ax, and as 7Ar is a constant, this leads to the decrease of 
the time step according to the relation At = '^j^{AxY / D. This property will be analysed 
in more detailed below. On the other hand, the increase of accuracy can be obtained, for 
fixed Ax and At, increasing the number of neighborhood cells in (2.6). 

3) If the cti are positive constants, we can interpret the connectivity coefficients ai in 
(2.6) as being proportional to the transition probability of a particle to jump from a cell to 
its neighborhood cells. More precisely, if Pi->j is the transition probability of a particle 
initially at [iAx — Ax/2, iAx + Ax/2] to be at [jAx — Ax/2,jAx + Ax/2], after a time 
At, by (2.6), we have 

N 

P,_>, = l-27ArV^ = l-^, P,_>, =7^^^L with|i-j|<Ar (2.22) 

if7iv<7iV, Tab. 2. 

The above numerical method can be easily extended for systems of nonlinear RD 
partial differential equations. To be more specific we take the RD equation 

f^Dg./(.) (2,23, 

11 



where f{(p) is a continuous function. The finite differences approximation to (2.23) can be 
written as 



N 

t-i-l h . 

V. 



fc+1 1- ^"^ "^ 



n + 7^ E i K- + <+^ - 2^') + ^^/(^') (2.24) 



i=l 



where, as before, Xl^ii = 1- By Theorem A, as •-jn is a scahng constant. At = '^n{Ax)'^/D 
and the error in the sup norm between the exact and approximate solutions of (2.23) and 
(2.24) is of the order of (Ax)'^. In appUcations, the finite difference approximation (2.24) is 
sufficient to obtain quahtative simulation results of RD equations. On the other hand, due 
to the characteristic finite scales in biological systems, some authors consider (2.24) as the 
basic RD equation [Turing, 1952]. However, in chemical systems we encounter stiff local 
kinetics. In some of this cases, the Euler type approximation (2.24) remains convergent. 

In order to compare the solutions of the diffusion and RD equations (2.1) and (2.23) 
with the solutions of the difference equations (2.6) and (2.24), we take the simplest in- 
tegrable reaction-diffusion system obtained by coupling the pseudo-first order reaction 
A + X -^^ 2X with diffusion, [Tilden, 1974], where X and A represent two chemical 
species and /3 is the reaction rate. Representing the chemical species and its concentra- 
tion by the same symbol, the RD equation for the pseudo-first order reaction in infinitely 

extended media is 

dX „ d^X 



at ^H^^"^'' f^'^^' 

where A and (3 are kept constants. We consider now the one dimensional lattice with 
coordinates (zAx). The initial distribution of X, localized at a; = 0, is represented by the 
function X{x, t = 0) = AxXq5{x), whose mean value in the interval [— Aa;/2, Aa;/2] is Xq. 
So, the solution of (2.25) with this initial distribution is 

X{x, t) = AxXq , e"" /^^* (2.26) 

We now approximate the reaction-diffusion equation (2.25) by the discrete evolution equa- 
tion 



N 

V. 



k+1 ^ ^"^ o; 



< + ^nJ2^ i<-^ + <+^ - 2^n) + ^t(3Av'n (2.27) 



n - , ^ 



where 7Ar = DAt/{Ax)'^ and the constants ai are given by Theorem A, Tab. 1. Under 
these conditions, 7Ar is a constant and At is determined by the relation At = 7Ar(Ax)^/i3>. 
The initial conditions corresponding to the localized distribution of X are then 

v^q=Xq, t;° = for z ^ 

Note that if X{x) is an initial distribution of matter and we take the class of func- 
tions yx = {y{x) : JiAx±Ax/2^(^)^^ = /iAx±Ax/2^(^)'^^'f°^ ^^^ ^ ^ ^}' *^^ numerical 
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solution of the diffusion equation at time t is the same for any Y G 3^x- Therefore, we 
define e = sup^f^ij^^^_^^/2^iAx+Ax/2] l^(^) " -^i\^ where Xi is the mean value of X{x) in 
the interval [iAx — Ax/2,iAx + Ax/2], and £ is z independent. The constant e is the 
incertitude in the initial function X{x) associated to the discretization step Ax. 

In Fig. 2 we show the global error, calculated with the norm of the supremum, 
between the solution (2.26) and the solution of (2.27), as a function of 7Ar, for A^ = 1 and 
A^ = 3. In the case A^ = 3, and in order to obey the relation ^ a^ = 1, the values of ct^, 
i = 1,2,3, are function of 73, obtained by solving (2.11a): ai = 3/2 — I373/4 + 673/2, 
a2 = -3/5 + 473 - 47I and 03 = 1/10 - 373/4 + 37I/2. 

We now estimate the computing time needed to calculate the solution of the discrete 
diffusion equation as a function of the number of neighbors taken into account into (2.6) 
and of the fixed time step At, for a finite lattice with zero fiux boundary conditions. If L 
is the number of integer sites of the lattice, the number of memory calls to calculate the 
concentration at time t is, disregarding boundary conditions, 

T^'^ - L{2N + l)Integer^ = L{2N + l)Integer — -^— (2.28) 

To decrease the global error four orders of magnitude, the increase in memory calls, or 
relative excess computing time, is 



^^+2 _ 2N + 3 7^ 



-,dif 



■n' 



With the values of Tab. 1, we have, T^j( = 0.79, T^jl = 0.84 and Ty^ = 0.88. Therefore, 
the computing time is faster with increasing N, with better accuracy in the global error. 
However, the incertitude e in the initial function X{x), introduced by the discretization 
step Ax, increases. This is a consequence of the scaling relating 77V =constant. However, 
for RD equations the global error can be larger for larger A^. In fact, the error associated 
to (2.24) is of the order of At = '^^{AxY / D and 7Ar+2 > 7Ar- 

To decrease the global error in the numerical integration of RD equations with in- 
creasing A^, we must therefore ask that, A-^+^t < A-^t, where A-^t is the time step taken 
in (2.24) with A'" neighbors. This last condition implies that 7Ar_|_2(A-^"'"^a:)^ < 7Ar(A^a:)^. 
With L(^+2)(A^+2x) = U^\A^x), where L^^) is the number of integer sites of the 
lattice for fixed A^, the relative excess computing time to decrease the global error four 
orders in the discrete RD equation is, by (2.28), 

.RD (L(^+2))3(2Ar + 3)7Ar ^ 2Ar + 3 /^^^ 



and 



^(^+^)/^" (LW)3(2Ar + 1)7^+2 - 2Ar + lV 7^ ^^"^"^ 



A(^+2)x < A(^)x J^^ (2.31) 

V lN+2 
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So, as 7Ar+2 > 7iv, the global error in RD systems decreases with increasing A^ if the 
lattice space increment decreases according to (2.31). With the values in Tab. 1, we have, 
A^^^x < 0.69A(^)x, A^^^x < O.SlA^^^x, A^^^x < O.SGA^^^a:. In this case, computing time 
is slower and we have, T^,^ = 2.41, T^,^ = 1.59 and Ty^ = 1.37. Therefore, in order to 
decrease the global numerical error in the integration of RD equations, it is sufficient to 
increase the number of neighborhood connectivities and simultaneously decrease the space 
lattice increment Ax according to the relation (2.31). 

For example, in the simulations of Fig. 2, the number of iterations k to reach the time 
t = 0.002 varies with 77V, k = Integer{t/ At) = Integer{Dt/'yN{Ax)'^). For A^ = 1 and 
A^ = 3 we have, k = 120 and k = 95, respectively. According to (2.29), the relative excess 
computing time is T^^ = 0.79, and therefore, Tg*-^ < T-^^^. For A^ = 3, to decrease the 
global error, we should have chosen Ax — 0.0069, according to (2.31). But, in this case, 
Tj^f = 2.41, and k = 289. 

2.2 - Two space dimensions 

To extend the previous results for the two-dimensional diffusion and RD equation, we 
consider a squared lattice with spatial coordinates (iAx^jAx) with z,j G Z. In order to 
enumerate lattice points in the neighborhood of the generic point {iAx,jAx), we associate 
to each integer n > 1 the set of integer coordinates J^ = {{r^s) : r^ + s^ = d^ ,r,s G 
Z , (i„ G N}, where {dn} is the sequence of positive integers such that r^ + s^ = d^ has 
integer solutions, {dn} = {1, 2, 4, 5, 8, 9, 10, 13, . . .}. For example, 

Ji = {(-1,0), (1,0), (0,-1), (0,1)} 

J2 = {(-1,-1), (1,-1), (-1,1), (1,1)} 

J3 = {(0,2), (2,0), (0,-2), (-2,0)} 

J4 = {(2, 1), (1, 2), (-2, 1), (-1, 2), (2, -1), (1, -2), (-2, -1), (-1, -2)} 

Let hn be the number of elements in each set Jn, hn = iJ^Jn- Writing the diffusion equation 
in two space dimensions as 

dct^ nv^ /a2</, a2</>\ 

where a^ are non negative constants with X]n=i <^n = I5 the finite difference approximation 

to (2.33) is 

M 

ar 



-■^r=vl+iiNY.Tir E «..+.- -y (2-34) 



n=l (r.,s)6J„ 



where the index n refers to the neighborhood order relative to the central cell with lattice 
coordinates (i,i), and 7Ar = DAt/{AxY, Fig. 3. The integer M is the space width of 
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the finite difference equation (2.34). In order to distinguish the error order index A^ and 
the neighbor order M we have introduced two different indices. The role of this indexes is 
specified in Theorem B below. 

The stabihty condition for (2.34) is easily derived (Sec. 2.1). The two-dimensional 
finite difference equation (2.34) is stable and obeys the maximum principle if 

and the constants a„ are non negative. 

The solutions of the finite difference equation (2.34) approach the solutions of the 
two-dimensional diffusion equations (2.33) under the following conditions: 

Theorem B. Let Ri be the rank of the matrix of the hnear system of inhomogeneous 
equations 

ml [k) (2my\2k ) 2^d^h. ^ / ' "" ^^'^^^ 

in the variables an, with, 1 <m < N + 1, < k < m/2, M = 1 + N + {N + 1)^/4, and 
let 7iv be as in Theorem A, for odd N. Let R2 be the rank of the matrix of the linear 
system of inhomogeneous equations (2.36), with M = Ri. Then, there exist real constants 
ai, . . . , ttM, solutions of (2.36), with M > R2 and R2 < 1 + A^ + (iV + 1) V4. Moreover, 
if M > R2, the linear system (2.36) has an infinite number of solutions. If, for some M, 
«!,..., aM can be chosen non negative and 7Ar < 7^, where 7^^ is given by (2.35), then 
the solution of the difference equation (2.34) approaches the solution of the continuous 
equation (2.33) with an error of the order of (Aa;)^*-'^"'"^-*. 

Remarks on Theorem B: It can happen that an error of the order of (Ax)"^^^^"^' 
occurs for an infinite set of solutions of (2.36), and to obtain the non negativity of the 
an, we must increase M. For example, if A^ = 3, i?i = 7 and M = 7, and On > for 
n = 1, . . . , 7. For A^ = 5, i?i = 14 and non negative solutions for the a^'s is obtained 
only if M = 15. In this case, it is possible to show that On > 0, for n = 1, . . . , 15 and 
ai5 G [3.7404 x 10~^, 4.747 x 10~^]. We have tested numerically the choice of positive 
constants ct^ up to A^ = 7. 

Proof: Let (/){x,y,t) be a solution of the two-dimensional diffusion equation (2.33). 
We suppose that (/)(a;, y, t) is analytic in x, y and t. Let us define the difference operator 

M 

ar 



A(l) = (l){x,y,t + At)-(p{x,y,t)-4'yN^-Tj- ^ {(p{x + rAx,y + sAx,t) - (l){x,y,t)) 

n=l "■ "■ (r,s)eJn 

(2.37) 
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where ^^=i ctn = 1 and •-jn = DAt/{Ax)'^. By the Taylor expansion, 

m=l k=0 ^ • \ / 

Now suppose that (r', s') G J^- By definition of J^ we have that (— r', s') G Jn and 
(r', — s') G Jn- So, with ei = rAx and 62 = sAx, if m is odd, Xl(r s)eJ„ r^~^s'^ = and 
we have 

y {(j){x + rAx.y + sAx.t) — (f){x,y,t)) 

{r,s)eJn. 

-\-oo 2m 



" ^ ^ (2m)! V A; J dx'^'^-^dy^ I ^ 

m=l fc=0 V ; \ / W \(r,s)eJ 



^2m-k^k 



{r,s)eJn 
2m— k „k 



Analogously, if k is odd, '^Zr seJ ^ '^ ~ ^' ^'^'^ 

^ ((/)(x + rAx, y + sAx, t) — (/)(a:, y, t)) 

(r,s)6J„ 

_ +2^ ^ (Aa:)2m /2^\ d^'^(j){x,y,t) I ^ 

m=\ fc=0 V ^ V / ^ \(r,s)6J 



(2.39) 



2m-2A: 2fc 



- 'J n 

Taking time derivatives of the two-dimensional diffusion equation, 

^ = D-fM ^"^"^ (240) 

introducing (2.39) and (2.40) into (2.37), and truncating the expansion for m = A^+ 1, we 
obtain for the operator A(f), 



N + l m ^2™, _! /m /„ X .. /..„ X M 



{r,s)eJn 
2(Af+2)' 



Ad^TiAx^'-T ^^ hirm\4j^r2m\y.^^ y. 

V Z^y -^1 A^ g 2m-2kg 2k \^i\f.i (2m)l\2kJ^ dnhn ^ 

m=\ k=0 \ \ / V ^ \ /n=l "■ " (r,s)e. 

+ O (a Ax 
Introducing the operator 



2m-2k 2k 



(2.41) 



M 



n=l (r,s)6J„ 



in order to make \A(j) — A^v\ = O ((Ax)^*-'^"'"^-'), we must have 

7^ /m\ ^ /2m\ ^ ^^ V r^— 2V^-0 (242) 
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withm = 1,...,A^+1 and A; = 0, ...,m. For fixed A^, we have ((A^ + 2)^ + (Ar + 2))/2 - 1 
equations and M unknown constants ai, . . . ,aM- However, the equations in (2.42) are 
not independent. Dependent equations in (2.42) occur for values of /c, say k' and k" , 
such that k' = m-k-, because (^) = (,-,). (^^) = &,) and E(.,.)eJ„ ^'"^-'''«''' = 
Sfr s)eJ r^"^~^'^ s^'^ . Therefore, we can restrict the range of the indix k in (2.42) to the 
integers k < m/2, obtaining (2.36). 

Let bm be the cardinahty of the set {{^) : < k < m/2}. By the Pascal triangle 
rule, it is straightforward to show by induction that, bm+2 = bm + '2, with bi = 1 and 
b2 = bs = 2. Therefore, for each m, the number of equations in (2.42), with k < m/2, is 
bm = 3/4 + (— 1)™'/4 + m/2. Now let us denote by Cm the number of equations in (2.42) 
up to order m and with k < m/2. So, we have the recurrence relation Cm,+i = c^, + ^m+i, 
with ci = 1. The solution of this recurrence is c^, = —1/8 + (— 1)"^/8 + to + m? /A, and 
the number of equations in (2.42), with k < m/2, is 

Ml - -^ + ^(-1)"^+' + iV + 1 + l(iV + 1)2 
o o 4 

Let Ri be the rank of the matrix of the linear inhomogeneous system (2.36), for A^ odd 
and M = Ml = A^ + 1 + (A^ + if /A. Let 7Ar be as in Theorem A. If i?i = Mi, then 
M = Ml = A^ + 1 + (A^ + 1)2/4, and the linear inhomogeneous system (2.36) has a 
solution. If, Ri < Ml, let R2 be the rank of the matrix of the linear inhomogeneous 
system (2.36), for A^ odd and M = Ri. Solving a linearly independent subsystem of (2.36) 
of dimension R2 x R2, for cti, . . . , an^, we obtain a solution that annulates R2 equations 
in (2.36). If we let M > R2 this linearly independent subsystem of (2.36) has an infinite 
number of solutions, dependent of M — R2 variables a^. But, in this case {Ri < Mi), 
these solutions also annulate the remaining dependent equations in (2.36). Otherwise, as 
the ctn's do not depend on 99, \A.ip\ = 0{Ax)'^^-^ ^'^\ with A^' < A^, and this contradicts 
Theorem A for the choice Lp{x, y, t) = Lp{x, t). 

Finally, if a^ > 0, for 1 < i < M and some M > R2, hy the maximum principle, 
the solution of the difference equation (2.34) approaches the solution of the continuous 
equation (2.33) with an error of the order of (Ax)2'-^+2) Therefore, Theorem B is proved. 

We show in Tab. 3 the values of the constants ai and 7Ar for A^ = 1 and A^ = 3 
(M = 3 and M = 7), for the discrete diffusion equation (2.34), as well as the values of 
stability limit (2.35). To obtain comparable results in dimensions 1 and 2, from the point 
of view of global errors, we must take the corresponding difference methods for the same 
value of A^. (Note that in dimension 1, M = A^). The scaling constant is the same for 
the same value of A^, but the connectivities with neighborhood sites must follow the order 
defined by Theorems A and B. 
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The connectivity relations determined by the above theorem have a simple geometric 
meaning. As a matter of fact, for each A^, the local connectivities with neighborhood 
sites, with a^ > in (2.34), follow a spherically symmetric front. Fig. 3, suggesting the 
possibility of elimination of lattice symmetries in numerical solutions of RD equations. 

To test the symmetry properties of the solutions of the discrete methods (2.34), we 
take as prototype model the Brusselator, (1.3). By (2.34), the difference RD equation for 
the Brusselator is 



M 
k A \ ^ '^^ \ ^ Ik k \ 



fe+1 

n=l {r,s)eJn 

+ At{kiA - {k2B + A;4X.,, + A;3«.,,) Vt,j) 
D ^ 

n=l l^r,s)eJn 



(2.43) 



+ AtihB^l^^^ - fc3(v'i.w)V2 



^,j^ ^^,^,j' 



where jn = ^maxAt/(Ax)^ and L'max = TXia,x{Di, D2}. 

In general, for systems of RD equations, the integration algorithm is a set of equations 
of type (2.43), one equation for each diffusive variable (pi. The consistency with the scaling 
relation 7^=constant, given by Theorem B, shows that we can only have D2 = Di or 
D2 = 0, in agreement, respectively, with the simulation and model choices of Keener & 
Tyson [1992] and Zhabotinsky & Zaikin [1973]. To maintain compatibility with these 
limits, the diffusion terms are weighted by the non dimensional factors -Di/-Dmax, where 
-Dmax = niax{Di, D2}. The space scale is now determined by Ax = v^maxAt/TJv, as 
otherwise we would be simulating diffusion below the mean free path of one of the diffusive 
variables. 

Under these conditions and for the same parameter values as in Fig. 1, we choose 
A^ = 1 (M = 3) in (2.43) and we have followed the time evolution of (/^i and (f2 with (2.43). 
In Fig. 4a) we represent the spatial distribution of Lpi at the time t = 435, for At = 1/6, 
in a square lattice of 525 x 525 sites, zero flux boundary conditions, and initial conditions 

ifl,^^=Aki/h, ipl,^^=Bk2k^/{Akiks), for (z,j) 7^ (262,262) 

(2.44) 

V?,262,262 = l-HAki/k^), V'2,262,262 = l-HBk2k4/ (Akiks)) 

The system develops undamped spherically symmetric travelling waves, originated at the 
perturbed lattice point {i,j) = (262,262). To analyse the radial symmetry of wave fronts, 
we measured the distance from a point in the wave front to the point where the perturbation 
has been initiated, as a function of the polar angle 9. We denote this distance by rg. In 
this case, the reference point at the wave front corresponds to a radial local maxima of 
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ipi- The normalized distance from the wave front to the center, rg/f, as a function of 
the polar angle 9, is represented in Fig. 5. From this numerical simulation we conclude 
that the maximum relative deviation from spherical symmetry of the wave fronts oscillates 
randomly in 9 with a maximum relative error of the order of 0.2%. The relative deviation 
as a function of 9 oscillates randomly, and the standard deviation of the fluctuations around 
the mean circular wave front is 1.2 x 10~^. Numerical analysis for linear RD equations 
show the same type of behavior for the global error as in Fig. 2. 

Other patterns with spiral symmetry that appear in experimental systems are plotted 
in Figs. 2.3b-d, [Agladze & Krinsky, 1982; Ross et a/., 1988]. 

2.3 - Three space dimensions 

For the three-dimensional diffusion equation, we consider a cubic lattice with spatial 
coordinates {iAx,jAx,kAx) with i,j,k G Z. We associate to each integer n > 1 the set 
of integer coordinates J„ = {(?, r, s) : g^ + r^ + s^ = (i„ , g, r, s G Z , (i„ G N}, where 
{dn} is the sequence of positive integers such that g^ + r^ + s^ = dn has integer solutions, 
{dr^} = {1, 2, 3, 4, 5, 6, 8, 9, 10, 11, . . .}. For example, 

Ji = {(-1,0,0), (0,-1,0), (0,0,-1), (0,0,1), (0,1,0), (1,0,0)} 
J2 = {(-1, -1, 0), (1, -1, 0), (-1, 1, 0), (1, 1, 0), (1, 0, 1), (-1, 0, 1), (2.45) 

(1, 0, -1), (-1, 0, -1), (0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1)} 

Writing the diffusion equation in three space dimensions as 

where J2n=i "^n = 1 and M is the space width of the connectivity constants ctj, the finite 
difference approximation to (2.46) is 

M 

\j:k = ^^,J,k + 67iV 2^ -fir ^ (^^+g,J+r,fc+s " ^^,J,fc) (2-47) 

n=l (q,r,s)eJn. 

where, hn = iJ^Jn, the index n refers to the order of neighborhood of each cell and 7Ar = 
DAt/(Aa;)^, Fig. 6. The three-dimensional finite difference equation (2.47) is stable and 
obeys the maximum principle if 

The solutions of the finite difference equation (2.47) approache the solutions of the 
three-dimensional diffusion equations (2.46) under the following conditions: 
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Theorem C. Let Ri be the rank of the matrix of the hnear system of inhomogeneous 
equations 

^mij \m2j (2m)! V2miy' \2m2) 4l d^hn , ^ , 

(2.49) 
in the variables an, with, l<m<N + l, 0< mi < m, < m,2 < nii, M = 3 + 
13N/3 + 3N^/2 + N^/6, and let 'Jn be as in Theorem A, for odd N. Let R2 be the 
rank of the matrix of the linear system of inhomogeneous equations (2.49), with M = Ri. 
Then, there exist real constants cti, . . . , aM, solutions of (2.49), with M > R2 and R2 < 
3 + 13A^/3 + 3A^2/2 + N^/6. Moreover, if M > R2, the linear system (2.49) has an infinite 
number of solutions. If, for some M, cti, . . . , aM can be chosen non negative and '-jn < Tjvj 
where 7]^ is given by (2.48), then the solution of the difference equation (2.47) approaches 
the solution of the continuous equation (2.46) with an error of the order of (Ax)^*^^"*"^-*. 

The proof of Theorem C is analogous to the one of Theorem B. 

In Tab. 4 we show the values of the constants a^'s, for A^ = 1 and A^ = 3. In the 
case of A^ = 3, the solution of (2.49) for positive ai can not be obtained for M = R2 = 9, 
and we have introduced the additional constant aiQ. So, system (2.49) has an infinite 
set solutions, as a function of aiQ. Numerical analysis shows that it is possible to choose 
ai > 0, with i = 1, . . . , 10, if ctio G [0.01278, 0.01785], and 73 < 73. In this case, the error 
between the solutions of (2.46) and (2.47) is of the order of (Ax)^'^, for any choive of ctio 
in the interval [0.01278,0.01785]. 

To test the symmetry properties of spherical wave fronts of the explicit method (2.47), 
for A^ = 1, we again consider the spatially extended Brusselator model with zero flux 
boundary conditions. We take as initial conditions a central perturbation to the steady 
state. In order to maintain the computations within the range of a personal computer, 
the three dimensional lattice has been chosen with physical dimensions of 175 x 175 x 175 
cells, and the time increment was At = 0.1. In Fig. 7 we present a three dimensional view 
of the concentration of the chemical species ipi. In this case, as we choose the time step 
as independent increment, the dimension of the simulation corresponds to a cube of side 
L = 175Ax = 175v/L'i At/71 ~ 136, where Di = 1.0 and D2 = 0. 

To analyse the symmetry properties of the numerically generated pattern in Fig. 7 and 
make the error analysis, we compare one-dimensional concentration profile extracted from 
one, two and three-dimensional simulations. If ip{x, t) is a solution of the one-dimensional 
RD equation, it is also a solution of the two and three-dimensional RD equations, with 
adapted initial conditions. Therefore, one-dimensional concentration profiles calculated 
numerically in the three cases should coincide. In Fig. 8 we present the three concentration 
profiles calculated with the finite difference methods (2.6), (2.34) and (2.47). 
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The matching of the concentration profiles shown in Fig. 8 show that the global 
errors between exact and numerical solutions of three-dimensional RD equations, as well 
as its symmetry properties, are the same as in the one and two dimensional cases analysed 
in Sec. 2.1 and Sec. 2.2. This fact enables to extrapolate properties of spherical and 
cylindrical solutions in two and three space dimensions from simpler and less computer 
time consuming one-dimensional simulations. 



3 - Conclusions 

We have presented a general class of finite difference approximations to diffusion and 
RD partial differential equations whose solutions approach the solution of the continuous 
system, within a global error of any prescribed order. We have shown that the differ- 
ence equations obey a scaling relation, relating space and time integration steps. This 
scaling relation is reminiscent from the well known invariance of the diffusion equation 
for transformations of the form x' — * ax and t' -^ a^t. The spurious lattice symmetries 
that appear in numerical simulations of RD equations have been controlled and minimized 
below reasonable errors. 

Numerical experiments. Fig. 1, show that specific patterns in RD systems are related 
with the temporal and spatial scales arising naturally in these systems. Therefore, any RD 
pattern obtained numerically without the proper calibration of the space and time scales 
can be qualitatively incorrect when compared with the exact solution of the continuous 
system of RD equations. 
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For computer calculation in single precision, the case A^ = 1 is in general sufficient to 
avoid lattice spurious symmetries. To be more specific, for the case A^ = 1 and global errors 
of the order of (Ax)^, the numerical alghorithms for the integration and RD equations are: 

dimension 1: 

^^' =^' + \i^ti + <i - 2^^f ) + At/(t^f ) 
dimension 2: 

dimension 3: (4.1aJ 

^I'^ + l ^ii'^ -I- (i)^ 4- ii'^ -U 1)^ -I- ii'^ -I- ii'^ 

i,j,m i,j,m '^ -i Q\^i—l,j,m '^ ^i-\-l,j,Tn '^ ^i,j — l,Tn '^ ^i,j-\-l,m '^ ^i,j,m — l 

'^ '^i,j,m+l ~ t'f^iJ,TnJ + TT7T\'^i-l,j-l,m + '^i+l,j-l,m + '^i-l,j + l,m 

+ ^i+l,j + l,m + ^i+l,j,m+l + "^i-l.j^m+l + '^i+l,j,m-l + ^i-l,j,m-l 
+ "^ij + lim+l + "^iJ-l.Tn+l "I" ^j,j + l,m-l + ^i,j-l,m-l ~ ^'^'^i,j,m) 

and the space and time steps are calibrated according the scaling relation 

^^* ' (4.16) 



(Aa;)2 6 

The finite difference methods (4.1) enable to measure relevant physical quantities in 
numerical simulations, in order to calibrate system parameters. The strategy of calibra- 
tion and validation of numerical methods leads to the possibility of analysing computer 
simulations as experiments with the same type of data analysis and conclusions. This is 
particularly important for modeling nonlinear RD systems, where analytical solution of 
transient processes are difficult or impossible to find. 

Several authors have developed different strategies to calibrate the numerical solutions 
of diffusion and RD equations. In particular, Dejak et al. [1987], based on the compari- 
son of numerical and known solution of diffusion equations found the heuristic condition 
7 = 1/6. In the context of cellular automata simulations Weymar et al. [1992] compared 
the preservation of symmetries of several Laplacian mask operators as a function of neigh- 
borhood connectivities. The results presented here give the precise foundation for these 
heuristic techniques. This is particularly important for the simulation and bench-marking 
of three-dimensional systems where complex topological structures are expected to appear. 
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On the other hand, in the two-dimensional case, the coefficients in (4.1a) are similar to 
those obtained in the numerical analysis of elliptic problems [Boisvert, 1981] which evolves 
the discretization of the Laplacian operator. 
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Table Captions 



Table 1: Parameters for the discrete diffusion equation (2.6), for A^ = 1,3 and 5, 
according to Theorem A. We have tested numerically the positiveness of the constants ai 
up to A^ = 23, for odd A^. In this case, for example, 023 — lO"^'^. Other values of •jn are: 
77 = 0.713841, 79 = 0.896295 and 711 = 1.07875. 

Table 2: Transition probabilities to neighborhood cells in the Brownian motion in- 
terpretation of diffusion. These transition probabilities depend on the scaling relation 
7Ar =constant and of the non negativity of the connectivity coefficients a^. 

Table 3: Parameters for the discrete diffusion equation (2.34), according to Theorem 
B, for A^ = 1 and A^ = 3. 

Table 4: Parameters for the discrete diffusion equation (2.47), for A^ = 1 and A^ = 3, 
under the conditions of Theorem C. In the case A^ = 3, equations (2.49) have an infinite set 
of solutions. The constants a^ with n = 1, . . . , 10, are positive if aio G [0.01278, 0.01785]. 
For A^ = 3, the parameters were obtained with the choice of ctio in the middle of the 
interval [0.01278,0.01785]. 
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Figure Captions 

Figure 1: Temporal evolution of non-equilibrium travelling waves obtained with the 
Brusselator model (1.3), for several values of the parameter 7 < 1/4. The initial condition 
has been chosen to produce spiral patterns at 7 = 0.01. The simulation parameters are 
Di = 0.008 , D2 = 0.0016, ki = 1.0, A = 1.0 and B = 2.3. We represent the concentration 
profiles for the chemical species (fi. Time evolution has been calculated on a lattice of 
200 X 200 cells, with zero flux boundary conditions, the same initial condition and a 
time step At = 0.1. To maintain the same radial velocity, the lattices have been scaled 
according to the relation Ax = ^J{AtDmaxh)■, where Dmax = max{il>i, D2}. Otherwise, 
propagation velocities would be different, and dependent on the integration space and time 
steps. 

Figure 2: Error between the solution (2.26) and the solution of (2.27), as a function 
of 7Ar, for A^ = 1 and A^ = 3, with the parameter values of Tab. 1. The error function 
is supj \X{iAxj kAt) — ff^|, calculated at the time t = 0.002. The simulations have been 
performed with (3 = 1.0, A = 1.0, D = 1, Ax = 0.01, At = -fN{Ax)'^/D and the initial 
value Xq = 0.01. For A^ = 1, the error is in fact minimized when 71 = 1/6 (Theorem 
A) and the global error is (0.01)^ = 10~^^. For A^ = 3, the global error is of the order 
of lO"^'^, which is below the limit of precision of computer calculation in single precision 
(-10-6). 

Figure 3: Space width and connectivity constants ct^ for the discrete version of the 
two dimensional diffusion equation (2.34). If follows from Theorem B that local connectiv- 
ities with neighborhood sites, with a^ > in (2.34), follow a spherically symmetric front, 
suggesting the possibility of elimination of lattice symmetries in numerical solutions of RD 
equations. 

Figure 4: a) Two dimensional patterns generated by the discrete difference RD 
equation (2.43) with A^ = 1, for the Brusselator, with initial conditions (2.44) and the 
parameters of Fig. 1. The depicted flgures have been calculated in a lattice of 525 x 525 
sites, with zero flux boundary conditions. The anisotropy analysis of the wave fronts in a) 
is presented in Fig. 5. From c) to d) we present other patterns generated with different 
initial conditions, qualitatively similar to patterns found in experimental systems. 

Figure 5: a) Anisotropy analysis of the pattern of Fig. 4a). Numerically obtained 
wave fronts deviate from spherical symmetry, within a maximum relative deviation of the 
order of 0.2%. The standard deviation of the fluctuations around the mean circular wave 
front is 1.2 x lO'^. 
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Figure 6: Space width M or neigiiborliood order for tiie discretization of the three 
dimensional diffusion equation (2.47), for A^ = 1 and A^ = 3. 

Figure 7: Three dimensional spherical symmetric travelling waves obtained the finite 
difference approximation (2.47) for the Brusselator model in three-dimensional extended 
media. The kinetic parameters in this simulation are the same as in Fig. 1, and the 
diffusion coefficients are Di = 1.0 and D2 = 0. 

Figure 8: Two and three dimensional plane wave solutions, and one-dimensional 
concentration profiles of the variable (fi, at t = 125, for the Brusselator model in extended 
media, calculated with one, two and three-dimensional RD equations. The parameters are 
the same as in Fig. 7, with A^ = 1. 
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